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ABSTRACT 

We present Submillimcter Array (SMA) observations of several deuterated species in the disk around 
the classical T Tauri star TW Hydrae at arcsecond scales, including detections of the DCN J=3-2 
and DCO + J=3-2 lines, and upper limits to the HDO 31,2—22,1, ortho-H2D + 1i,q— li i and para-D2H + 
li,o - lo,i transitions. We also present observations of the HCN J=3-2, HCO + J=3-2 and H 13 CO + 
J=4-3 lines for comparison with their deuterated isotopologues. We constrain the radial and vertical 
distributions of various species in the disk by fitting the data using a model where the molecular 
emission from an irradiated accretion disk is sampled with a 2D Monte Carlo radiative transfer code. 
We find that the distribution of DCO + differs markedly from that of HCO + . The D/H ratios inferred 
change by at least one order of magnitude (0.01 to 0.1) for radii <30 AU to >70 AU and there is a 
rapid falloff of the abundance of DCO + at radii larger than 90 AU. Using a simple analytical chemical 
model, we constrain the degree of ionization, x(e-)=n(e-)/n(H2), to be ~ 10 -7 in the disk layer(s) 
where these molecules are present. Provided the distribution of DCN follows that of HCN, the ratio of 
DCN to HCN is determined to be 1.7±0.5 x 10 -2 ; however, this ratio is very sensitive to the poorly 
constrained vertical distribution of HCN. The resolved radial distribution of DCO + indicates that in 
situ deuterium fractionation remains active within the TW Hydrae disk and must be considered in 
the molecular evolution of circumstellar accretion disks. 

Subject headings: circumstellar matter — comets: general — ISM: molecules — planetary systems: pro- 
toplanetary disks — stars: individual (TW Hydrae) — stars: pre-main-sequence 



1. INTRODUCTION 

Millimeter-wave interferometers have imaged the gas 
and dust surrounding over a dozen T Tauri and Herbig 
Ae stars (see Dutrey et al. 2007 for a review). These 
studies demonstrate the potential to dramatically im- 
prove our understanding of disk physical and chemical 
structure, providing insights that will ultimately enable 
a more comprehensive understanding of star and planet 
formation. As analogs to the Solar Nebula, these circum- 
stellar disks offer a unique opportunity to study the con- 
ditions during the planet formation process, especially 
the complex chemical evolution that must occur. In the 
outer parts of the disk directly accessible to millimeter- 
wave interferometry, observations of deuterated species 
are particularly important because they can constrain 
the origin of primitive solar system bodies such as comets 
and other icy planetesimals. 

Deuterated molecule chemistry is sensitive to the tem- 
perature history of interstellar and circumstellar g as, as 
well as to density-sensitive processes such as molecular 
depletion in the cold (< 20 K) and dense (> 10 6 cm~ 3 ) 
disk midplane. The similarity of molecular D/H ratios 
between comets and high mass hot cores has been used 
to argue for an "interstellar" origin of cometary mat- 
ter, but ambiguity remains and this argument is not se- 
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cure (Bergin et al. 2007). Aikawa & Herbst (1999) in- 
vestigate the chemistry of deuterium-bearing molecules 
in outer regions of protoplanetary disks and find that 
molecules formed in the disk have similar D/H ratios 
as those in comets. To date, the determination of D/H 
ratios in disks has been limited to the measurements of 
DCO+/HCO+ in two classical T Tauri Stars (cTTs) with 
single dish telescopes (TW Hydrae, van Dishoeck et al. 
2003; DM Tauri, Guilloteau et al. 2006), but DCO+ has 
not been observed in comets. Although H2D 4 " and HDO 
lines have been detected in disks (Ceccarelli et al. 2004, 
2005), there are no corresponding H^~ and H2O obser- 
vations that allow derivations of the D/H ratio. There- 
fore, direct measurements in disks of species observed in 
comets, such as DCN and HCN, will be important to help 
relate deuterium fractionation in disks to that in comets. 

Deuterated molecules in cold pre-stellar and protostel- 
lar cores are found to be enhanced by orders of magnitude 
over the elemental D/H abundance ratio of 1.5xl0 -5 
through fractionation in the gas phase at low tempera- 
tures, an effect driven primarily by deuterated Hjj", which 
is formed by exchange reaction between Hjj~ and HD, 
then transfers its deuteron to neutral species. Theoret- 
ical models of disks (Aikawa et al. 2002; Willacy 2007) 
with realistic temperature and density structure show 
that the DCO + /HCO + ratio increases with radius due 
to the decreasing temperature moving out from the star. 
Spatially resolved data of deuterated molecules will help 
test disk physical models through these chemical conse- 
quences. 

The current capabilities of millimeter-wave observato- 
ries are limited both by sensitivity and by the small an- 
gular size of circumstellar disks, which makes spatially 
resolving the deuterium fractionation difficult. In this 
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work we take advantage of the proximity of the most 
nearby classical T Tauri star, TW Hydrae, which is sur- 
rounded by a disk of radius 3'.'5, or 200 AU at a distance 
of 56 pc (Qi et al. 2004), and a possible orbiting planet 
of mass (9.8±3.3) Mj upit er at 0.04 AU (Setiawan et al. 
2008), to study the physical and chemical structure of a 
protoplanetary environment at high spatial and spectral 
resolution using intcrfcromctry. The TW Hydrae disk is 
viewed nearly face-on, and so is well posed to investigate 
the radial distribution of molecular emission. Although 
the current angular resolution of <~l-2" places only a 
few independent beams ( or 'pixels') across the disk, ef- 
fectively improved resolution of the disk chemistry can 
be achieved thanks to the fact that the gas kinemat- 
ics are essentially Keplerian. Beckwith & Sargent (1993) 
show that molecular line emission from Keplerian disks 
displays a "dipole" pattern near the systematic veloc- 
ity due to the shear created by the orbital motion. The 
emission near the systemic velocity is dominated by the 
outer disk regions, and the separation of the emission 
peaks depends sensitively on the abundance gradient of 
the emitting species with radius. This important feature, 
commonly seen in velocity channel maps of disks with a 
range of inclination angles (even as small as 6-7 degrees 
in the case of TW Hydrae, Qi et al. 2004), may be used 
to study the radial distribution of molecular emission far 
more precisely than is afforded by the limited number of 
"pixels" provided by the available resolution. 

Here we report on Submillimcter Array (SMA) 5 (Ho 
et al. 2004) observations of deuterated species in the disk 
around TW Hydrae, including the first detection and im- 
ages of DCN and DCO + . In §2 we describe the observa- 
tions, while in §3 we introduce the analysis method used 
and the molecular distribution parameters derived. We 
describe the model fitting results and discuss the impli- 
cations in §4, and we present a summary and conclusions 
in §5. 

2. OBSERVATIONS 

All of the observations of TW Hydrae (R.A.: 
ll /l 01 m 51. s 875; Dec: -34°42'17."155; J2000.0) were 
made between 2005 February and 2006 December using 
the SMA 8 antenna interferometer located atop Mauna 
Kea, Hawaii. Table 1 and Table 2 summarize the ob- 
servational parameters for the detected and undetected 
species, respectively. The SMA receivers operate in a 
double-sideband (DSB) mode with an intermediate fre- 
quency band of 4-6 GHz which is sent over fiber op- 
tic transmission lines to 24 "overlapping" digital correla- 
tor "chunks" covering a 2 GHz spectral window. Two 
settings were used for observing the HCN/DCN and 
HCO+/DCO+ lines: at 265.7-267.7 GHz (lower side- 
band) the tuning was centered on the HCN J=3-2 line at 
265.8862 GHz in chunk S22 while the HCO+ J=3 -2 line 
at 267.5576 GHz was simultaneously observed in chunk 
S02; at 215.4-217.4 GHz (lower sideband) the tuning was 
centered on the DCO+ J=3-2 line at 216.1126 GHz in 
chunk S16 while the DCN J=3-2 line at 217.2386 GHz 
was simultaneously observed in chunk S02 (the HDO 
31,2-22,1 line at 225.90 GHz was also covered in the USB 

5 The Submillimeter Array is a joint project between the Smith- 
sonian Astrophysical Observatory and the Acadcmia Sinica Insti- 
tute of Astronomy and Astrophysics, and is funded by the Smith- 
sonian Institution and the Acadcmia Sinica. 



in chunk S06). Combinations of two array configurations 
(compact and extended) were used to obtain projected 
baselines ranging from 6 to 180 meters. Cryogenic SIS 
receivers on each antenna produced system temperatures 
(DSB) of 200-1400 K at 260 GHz and 100-200 K at 210 
GHz. Observations of ortho-H 2 D+ and H 13 CO + were 
simultaneously carried out (in the dual receiver observa- 
tional mode) on 28 December 2006 using only the com- 
pact array configuration. Three out of eight antennas 
were equipped with 400 GHz receivers (DSB system tem- 
perature between 400 and 800 K) and were available for 
observations of the 372.4213 GHz ortho-H 2 D+ li, -li,i 
line, and all eight antennas were operating with 345 GHz 
receivers (DSB system temperature between 150 and 300 
K) for observations of the 346.9985 GHz H 13 CO+ J=4-3 
and 345.796 GHz CO J=3-2 lines. The correlator was 
configured with a narrow band of 512 channels over the 
104 MHz chunk, which provided 0.2 MHz frequency res- 
olution, or 0.28 km s" 1 at 217 GHz, 0.23 km s" 1 at 267 
GHz and 0.16 km s _1 at 372 GHz. Observations of the 
691.6604 GHz para-D 2 H+ li,o-lo,i nne were shared with 
observations of the CO J=6-5 line made on 17 February 
2005 (details provided in Qi et al. 2006). Calibrations 
of the visibility phases and amplitudes were achieved 
with interleaved observations of the quasars J1037-295 
and Jl 147-382, typically at intervals of 20-30 minutes. 
Observations of Callisto provided the absolute scale for 
the calibration of flux densities. The uncertainties in the 
flux scale are estimated to be 15%. All of the calibra- 
tions were done using the MIR software package 6 , while 
continuum and spectral line images were generated and 
CLEANed using MIRIAD. 

3. SPECTRAL LINE MODELING 

Several model approaches have been developed to in- 
terpret interferometric molecular line and continuum ob- 
servations from disks (Dutrey et al. 2007). In brief, the 
approach of Qi et al. (2004, 2006) works as follows: the 
kinetic temperature and density structure of the disk is 
determined by modeling the spectral energy distribution 
(SED) assuming well mixed gas and dust, with the re- 
sults confirmed by the resolved (sub)mm continuum im- 
ages. Then, a grid of models with a range of disk param- 
eters including the outer radius R out , the disk inclina- 
tion i, position angle P.A. and the turbulent line width 
vturb are produced and a 2D accelerated Monte Carlo 
model (Hogerheijde & van der Tak 2000) is used to cal- 
culate the radiative transfer and molecular excitation. 
The collisional rates are taken from the Leiden Atomic 
and Molecular Database 7 (Schoier et al. 2005) for non- 
LTE line radiative transfer calculations. Specifically in 
the models used in this paper, the rate coefficients for 
HCO+ and DCO + in collisions with H 2 have been cal- 
culated by Flower (1999); the rate coefficients for HCN 
and DCN in collisions with H 2 have been scaled from the 
rate of HCN-He calculated by Green & Thaddeus (1974). 
The model parameters are fitted using a x 2 analysis in 
the (u, v) plane. 

In the studies by Qi et al. (2004, 2006), the distribution 
of CO molecules was assumed to follow the H nuclei as 
derived from the dust density structure and a gas of solar 

6 http://www.cfa.harvard.edu/~cqi/mircook.html 

7 http://www.strw.leidenuniv.nl/~moldata 
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TABLE 1 

Observational Parameters for SMA TW Hydrae (detected species) 

HCN 3-2 HCO+ 3-2 H 13 CO+ 4- -3 DCN 3-2 c DCO+ 3-2 

Rest Frequency (GHz): 265.886 267.558 346.999 217.239 216.113 

Observations 2005 Mar 04 2005 Mar 04 2006 Dec 28 2006 Apr 28 2006 Apr 28 

2005 Apr 21 2005 Apr 21 2006 Feb 03 

2005 April 26 2005 Apr 26 

Antennas used 7 7 8 8 8 

Synthesized beam: 1'.'6 x l'.'l PA -0.5° 1'.'6 x l'.'l PA -6.3° 4'.'1 x l!'8 PA 3.3° 5'.'9 x 3'.'2 PA -1.5° 2'.'6 x 1'.'6 PA 2.8° 

Channel spacing (km s" 1 ): 0.23 0.23 km 0.70 0.56 0.28 

R.M.S. a (Jy beam- 1 ): 0.35 0.29 0.16 0.10 0.10 

Peak intensity'' (Jy) 47 2M O70 031 056 

a SNR limited by the dynamic range. 

b Intensity averaged over the corresponding beam. 

c Only compact configuration data used for DCN observation. 



TABLE 2 

Observational Parameters for SMA TW Hydrae (undetected species) 



HDO 3i, 2-2 2 ,i o-H 2 D+ li,o-li,i p-D 2 H+ li,o-lo,i 



Rest Frequency (GHz): 
Observations: 
Antenna used: 
Synthesized beam: 
Channel spacing (km s -1 ): 
R.M.S. (Jy beam" 1 ): 

^molecule (1 g) (cm~ 2 ): 



225.897 
2006 April 28 
8 

5^' 7 x 3'.'1 PA -0.6° 
0.54 
0.11 
<2.0 x 10 14 



372.421 
2006 Dec 28 
3 

4'.'7 x 3'.'8 PA 14.5° 
0.65 
1.17 
<1.7 x 10 12 



346.999 
2005 Feb 17 
4 

3'.'3 x l'/3 PA 7.5° 
0.35 
7.39 
<9.0 x 10 14 



composition. However, molecules in disks do not neces- 
sarily share the same distribution as molecular hydrogen. 
Theoretical models (e.g. Aikawa et al. 1996; Aikawa & 
Nomura 2006) predict so-called three-layered structure; 
most molecules are photo-dissociated in the surface layer 
of the disk and frozen out in the mid-plane where most 
of hydrogen resides, with an abundance that peaks in the 
warm molecular layer at intermediate scale heights. To 
approximate this more complex behavior, we introduce 
new molecular distribution parameters for use in spectral 
line modeling. 

For first-order analysis in the radial molecular distri- 
butions, the column densities are assumed to vary as a 
power law as a function of radius: 

Ei(r) = £j(10AU)( 



10AU' 



Here £j and Pi describe column density distribution of 
a specific species (i) rather than disk (hydrogen) column 
density. In most disk models to date, a single radial 
power-law is adopted for fitting the hydrogen or dust sur- 
face density. By assuming CO follows the distribution of 
hydrogen column density, Qi et al. (2004, 2006) also find 
satisfactory power-law fits for multiple CO transtions. 
When a single power law is insufficient for fitting the 
radial distribution, then a broken power law with two 
different indices will be used. 

To calculate the surface density at different radii, the 
vertical molecular distribution is needed, but this dis- 
tribution may well vary with distance from the star. 
However, theoretical models (Aikawa & Nomura 2006) 
indicate that the vertical distribution of molecules at 
different radii is similar as a function of the hydrogen 
column density measured from the disk surface, £21 = 
Eh/(1-59 x 10 21 cm -2 ), where the denominator is the 
conversion factor of the hydrogen column density to A v 
for the case of interstellar dust. As indicated by Figure 




0.001 0.010 0.100 1.000 10.000 100.000 

Fig. 1. — This schematic diagram shows an arbitrary molecular 
vertical distribution as a function of E21 measured from the disk 
surface at a certain radius. a m and a s are the midplane and surface 
boundary parameters for model fitting (see §3 of the text). 

5 of Aikawa & Nomura (2006), the vertical distribution 
of molecular abundances shows a good correlation with 
£21- Figure 1 shows a schematic diagram of an arbi- 
trary distribution of a molecule in disk where the x-axis 
shows £21 and the y-axis shows the normalized molecu- 
lar fraction. Smaller £21 values (to the left of the plot) 
approach the surface of the disk, while larger £21 values 
(to the right) approach the midplane of the disk. For 
modeling, we make the further simplifying assumption 
that gaseous molecules exist with a constant abundance 
in layers between a s and a m , the surface and midplane 
boundaries of £21, respectively, as indicated in Figure 1 
by the shaded area. While the vertical distribution of 
molecules in disks may have a more complex distribu- 
tion than assumed here, the adopted parameters provide 
a gross approximation to the vertical location where the 
species is most abundant. This is adequate for a first 
description, given the quality of the data available. For 
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Fig. 2.— Top: The HCO+, DCO+, HCN and DCN J=3-2 spectra at the peak continuum (stellar) position. The fluxes of HCO+ and 
DCO+ are averaged over the beam of DCO+ J=3-2 (2'.'6 X 1'.'6 PA 2.8°). The fluxes of HCN and DCN are averaged over the beam of 
DCN J=3— 2 (5'.'9 X 3"2 PA -1.5°). The vertical dotted lines indicate the positions of the fitted Vlsr f° r each molecular transition except 
for DCN 3-2 where we adopt the V LS R from that of HCN J=3-2. Bottom: Velocity channel maps of the HCO+, DCO+, HCN and DCN 
J=3-2 emissions toward TW Hydrae. The angular resolutions arc 1'.'6 X l'.'l at PA -6.3° for HCO+ J=3-2 and 1'.'6 X l'.'l at PA -0.5° for 
HCN J=3— 2. The cross indicates the continuum (stellar) position. The axes are offsets from the pointing center in arcseconds. The fa- 
contour steps arc 0.4, 0.12, 0.35 and 0.09 Jy beam -1 for HCO+, DCO+, HCN and DCN J=3-2 respectively and the contours start at 2a. 



example, the model of Aikawa & Nomura (2006) shows 
two vertical peaks of HCO + , but the secondary peak is 
an order of magnitude smaller, and provides effectively 
negligible emission. Also, at least for the nearly face-on 
disk of TW Hydrae, the uncertainties in vertical distri- 
butions do not affect the derived radial distribution. 

This simple model captures the basic characteristics of 
three-layered structure predicted by theoretical models. 
Using the new distribution parameters: Ej(10AU) and 
Pi (radial), <r s and a m (vertical), the radial and vertical 
distributions of the molecules in disks can be constrained 
by observation. The best fit model is obtained by mini- 



mizing the x 2 , the weighted difference between the real 
and imaginary part of the complex visibility measured at 
the selected points of the (u, v) plane. The x 2 values are 
computed by the simultaneous fitting of channels cover- 
ing LSR velocities from 1 to 4 km s~ . R out and pi are 
calculated on grids in steps of 5 (AU) and 0.2, respec- 
tively. log(<7 s ) and log (<r m ) are calculated on grids in 
steps of 0.2 within the range from —2 to 2. Xj(10AU) and 
i are found with each pair of radial and vertical distribu- 
tion parameters by minimizing x 2 ■ The systemic velocity 
Vls-R lS fit separately since it is not correlated with those 
distribution parameters. For each fit parameter, the la 
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uncertainties are estimated as Xia 



Xr, 



f2n, where 



n is the number of degrees of freedom and \m i s the % 2 
value of best-fit model, as in Isella et al. (2007). We 
tested values for the turbulent velocity in the range from 
0.0 to 0.15 km s _1 and found that exact value does not 
have a significant impact, in part because of the coarse 
spectral resolution of the data. Therefore, we fixed the 
turbulent velocity at an intermediate value, 0.08 km s" 1 . 



4. RESULTS AND DISCUSSION 
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Fig. 3. — The solid and dotted contours show the temperature 
and density profiles from the TW Hya model of Calvet et al. (2002). 
The blue and red lines confine the locations of HCO+ and HCN in 
the model (see text). 

Figure 2 shows the spectra of HCO+, DCO+, HCN 
and DCN at the stellar (peak continuum) position and 
their velocity channel maps. Table 3 summarizes the 
power-law fitting results on HCO+, DCO+ and HCN. 
The DCN and H 13 CO + lines are too weak for a \ 2 anal- 
ysis of their distribution parameters, and so only the 
ratios of DCN/HCN and HCO+/H 13 CO+ are fit. Fig- 
ure 3 shows the density and temperature contours of the 
disk model (adopted from Calvet et al. 2002) and the 
locations of HCO + and HCN derived from the model fit- 
ting procedure. Figure 4 shows the radial distribution of 
the molecular column densities of the best fit models for 
HCO+, DCO+ and HCN. Not surprisingly, we found that 
the radial distributions are better constrained than the 
vertical ones: the local minimum of \ 2 f° r different grids 
of vertical parameters are within the noise limit, i.e. the 
vertical parameters are the least well determined, due in 
part to the face-on nature of the TW Hydrae disk. We 
therefore treat the best-fit vertical results as fixed, and 
investigate the uncertainty of other parameters. 

4.1. HCO+ and DCO+ 

The first detection of DCO + in the disk around 
TW Hydrae was obtained by van Dishoeck et al. (2003) 
with the James Clerk Maxwell Telescope (JCMT), who 
reported a beam-averaged DCO + /HCO + abundance ra- 
tio of 0.035. Our spatially resolved observations of the 
DCO + and HCO + J=3-2 emission from the disk suggest 
a more complex chemical situation. 

Figure 5 shows the channel maps of HCO + J=3-2, to- 
gether with the best-fit model and the data-model resid- 
uals. Table 3 lists the best-fit model parameters. The \ 2 



surface for the radial power index Phco+ an( A the outer 
radius R OM t is shown in the top panel of Figure 6. The 
lcr contour confines R ou t to 200±10 AU and Phco+ to 
-2.9±0.3. The H 13 CO+ J=4-3 line was also detected, 
but the emission is not strong enough to constrain its dis- 
tribution. Assuming that H 13 CO + has the same distribu- 
tion as HCO+, then fitting the ratio of HCO+/H 13 CO+ 
to match the intensity of the H 13 CO + emission indicates 
an HCO+ /H 13 CO+ ratio of 100±25. This value is consis- 
tent with the nominal solar system value of 89. Figure 7 
presents the best-fit model spectra of H 13 CO + J=4-3 
overlayed on the SMA data. 

The vertical distribution of DCO + is even less well con- 
strained than that of HCO + , and we choose to adopt the 
same values of a s and a m for both HCO + and DCO + 
(as indicated in Table 3). Because of the nearly face-on 
viewing geometry for the disk, the fitting of radial distri- 
bution power-law index p^ is not affected significantly by 
the uncertainties in vertical structure, but only the value 
of £j(10AU) needs to be adjusted. 

Examination of the channel maps shows upon inspec- 
tion that the radial distributions of DCO + and HCO + 
are different. To demonstrate how differences in the ra- 
dial distribution affect the resulting images, Figure 8 
presents three models of DCO + radial distribution and 
Figure 9 shows the corresponding simulated channel 
maps of the DCO+ J=3-2 line. 

Model 1 assumes that the DCO + distribution fol- 
lows the best- fit model of HCO + . The minimum 
X 2 is determined to be 575396 and the corresponding 
DCO+/HCO+ is found to be 4.7 x 10" 2 . Comparing the 
simulated maps from Model 1 with the data in Figure 9 
shows a distinct difference in that the double peaked na- 
ture of the central channel in the data is more obvious 
than in this model. The contrast of the contour levels be- 
tween the central channel and the adjacent channels are 
also smaller in the data than in this model. Because the 
emission of the central channel mostly originates at large 
disk radii, these differences suggest that the DCO + emis- 
sion arising from the outer regions of the disk is stronger 
than predicted by this model. The slope of radial DCO + 
distribution does not decrease as steeply as does HCO + . 
In other words, the D/H ratio must increase with increas- 
ing radius. 

Model 2 shows the best-fit result for the radial dis- 
tribution of DCO + assuming a single power index. As 
shown in Figure 4, the radial distribution of DCO + is 
strikingly different from that of HCO + , with a positive 
power index of 2.4 and a smaller but better constrained 
Rout 90±5 AU. The contours of the iso-% 2 surface for 
DCO + in the middle panel of Figure 6 indicate that the 
uncertainty of the radial power index is very large but 
that the index is still larger than 1.6 within the la er- 
ror, much larger than —2.9 found for HCO + . The simu- 
lated DCO + channel maps for Model 2 shown in Figure 9 
are an improvement over Model 1 in matching the data. 
However in this model, R ou t (~ 90 AU) is much smaller 
than the disk radius observed with HCO + and CO (~ 
200 AU) and DCO + increases with radius but then disap- 
pears sharply at R ou t^90 AU, as a step function, which 
is hard to understand. Comparison with the data shows 
that there are fewer complete contours around the dou- 
ble peaks in the central channel, which suggests that the 
DCO + emission is maximum at an intermediate radius 
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Fig. 4. — Radial distribution of molecular column densities and DCO + /HCO + ratio for the best-fit models for TW Hydrae. Solid lines 
depict single power-law fits, while the dashed lines are for DCO + Model 3 where two power-law indices are used. 



TABLE 3 
Fitting results 



Parameters 


HCO+ 


DCO+ a 


DCO+ 6 


HCN 


Stellar Mass: M*(M ) 


0.6 


0.6 


0.6 


0.6 


Inclination: i(deg) 


6.8±0.3 


7.4±0.6 


7.4 


6.6±0.8 


Systemic velocity: V^sfl(km s —1 ) 


2.88±0.05 


2.94±0.06 


2.9 4 


2.73±0.06 


Position angle: P.A.(deg) 


-27.4 


-27.4 


-27.4 


-27.4 


Turbulent line width: V[„ r j(km s" 1 ) 


0.08 


0.08 


0.08 


0.08 


Outer radius: R out (AU) 


200 ± 10 


90 ± 5 


160 


100 ± 10 


Column Density at 10AU: S(10AU) (cm" 2 ) 


3.8 ±0.5 x 10 15 


1.9±0.2 x 10 10 


4.8 x 10 6 


2.4±0.4 x 10 14 


Radial power index: p 


-2.9 ± 0.3 


2.4 ± 0.8 


7,-6 c 


-1.0 ± 1.2 


Vertical Parameters: <r s ,<Tm 


0.1, 10 


0.1, 10 


0.1,10 


0.3, 30 


Minimum x 2 '- 


329588 


575347 


575347 


297198 


Reduced x 2 '- Xr 


1.21 


1.95 


1.95 


1.35 



a DCO+ Model 2. 

h DCO + Model 3, no error estimation. 
c DCO+ turning point at radius 70 AU. 

rather than at the edge. 

For Model 3, instead of using a single power index (p) 
to fit for the radial distribution of DCO + column density 
N(DCO + ), the fit uses two power indices (pi and p2) and 
a turning point (Rt) where the power law index changes 
from pi to p2, i.e. the location of the peak of N(DCO + ). 
The parameters pi, p2 and Rt are searched within lim- 
ited grids to minimize x 2 . In this model, N(DCO + ) is 
found to increase with radius out to ~ 70 AU, and then 
to decrease. Table 3 presents the best-fit parameters; 
no error estimation is provided due to the computation 
difficulties. Figure 9 shows the best-fit model images. 
Comparison of Model 3 with the data shows a visual im- 
provement over Model 2 in the central channel, although 
the x 2 value of the two models are not distinguishable. 
We thus cannot clearly discriminate with the \ 2 statistic 
if there is indeed a peak with radius for N(DCO + ) as 



in Model 3, or if DCO + increases with radius and dis- 
appears suddenly around 90 AU as in Model 2. Even 
with this ambiguity however, both models imply that 
the D/H ratios change by at least an order of magni- 
tude (0.01 to 0.1) from radii <30 AU to >70 AU and 
that there is a rapid falloff of N(DCO + ) at radii larger 
than 90 AU. Because the emission in the central channel 
comes from the outer part of the disk (Keplerian rota- 
tion) projected along the line-of-sight with a very small 
inclination of around 7 degrees for TW Hydrae, the cen- 
tral velocity channel is most important for constraining 
the radial distribution. Based on the difference in the 
central velocity channel map between the models, we be- 
lieve Model 3 is the more plausible description of the 
radial distribution of DCO + . Figure 10 shows the chan- 
nel maps of the DCO + J=3-2 emission, together with 
those of Model 3 and the data-model residuals. 
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Fig. 5. — Top: Velocity channel maps of the HCO + J=3— 2 emission toward TW Hydrae. The angular resolution is 1"6 X l'.'l at PA 
-6.3°. The cross indicates the continuum (stellar) position. The axes are offsets from the pointing center in arcseconds. The la contour 
step is 0.4 Jy beam -1 and the contours start at 2<r. Middle: channel map of the best-fit model with the same contour levels. Bottom: 
difference between the best-fit model and data on the same contour scale. 



Observations of deuterated molecular ions and the level 
of deuterium fractionation have been used to estimate the 
ionization degree in molecular clouds, and a similar anal- 
ysis can be applied to circumstellar disks. If we consider 
only the ionization balance determined by HCO + , H^, 
DCO + and electrons in steady state as shown in Equa- 
tion 14 of Caselli (2002), the electron fractional abun- 
dance can be derived to be around 10~ 7 . Of course, 
this value is only valid in the intermediate layer where 
HCO + is abundant and multiply deuterated Hjj~ is less 
abundant than HCO + . Several important complications 
are also neglected in this analysis, including the presence 
of other atomic and molecular ions, neutral species be- 
sides CO which destroy H2D+, and negatively charged 
dust grains and refractory metals. Still, accurate mea- 
surements of DCO + and HCO + are the first steps toward 
an understanding of the ionization fraction in the disk. 

The increase of D/H ratios from inner to outer disk is 
generally consistent with the current theoretical models 
of the gas-deuterium fractionation processes that con- 
sider the effect of cold temperature. But the quick dis- 
appearance of DCO + at radii beyond 90 AU (comparing 
with R ou t around 200 AU for CO and HCO + ) is puzzling, 
since DCO + is expected to be abundant and observable 
in the cold outer region of the disk where HCO + is still 
available. More theoretical work is needed to explain the 
disappearance of DCO + in the outer part of the disk. 

4.2. HCN and DCN 

Figure 11 shows the HCN J =3-2 channel maps, to- 
gether with the best fit model and residuals. Table 3 
lists the best-fit model parameters, and Figure 4 shows 
the radial distribution of column density derived. The 



X 2 surface shown in the bottom panel of Figure 6 indi- 
cates R out 100±10 AU and p HC N -1.0±1.2. The radial 
power index is poorly constrained probably due to more 
complex distributions for HCN. A detailed comparison 
of the molecular distributions of HCN and CN will be 
presented elsewhere. 

Although the best fit vertical parameters seem to in- 
dicate that HCN (ff m .HCN = 30) is found much deeper 
toward the midplane than is HCO + (a m .hco+ = 10)) we 
emphasize that we are not able to accurately constrain 
the vertical distributions from the present data. This am- 
biguity strongly affects the column density of HCN (i.e. 
Shcn(IOAU)) needed to fit the data (not the power in- 
dex phcn of radial distribution). A worse fit to the data 
(x 2 larger by 3c over the best-fit model) can be obtained 
by assuming that the vertical distributions of HCN and 
HCO+ are the same, but the HCN column density is 1.5 
times larger than that needed for the best fit model due 
to higher density near the midplane. 

The DCN J=3-2 transition is detected at a signal-to- 
noise ratio of 3 near the fitted HCN Vlsr of 2.73 km s _1 
(Figure 2). While this signal-to- noise ratio is not high, 
the significance of the detection is further supported by 
the channel maps (Figure 12 upper panel) where the ve- 
locity gradient along the disk position angle of ~ —30° 
is consistent with that seen in CO J=2-l and J=3-2 (Qi 
et al. 2004) and the other molecular lines presented in 
this paper. Since the DCN 3-2 emission is weak, we 
are not able to fit for the molecular distribution and so 
make the simplifying assumption that the distribution of 
DCN follows that of HCN. As with H 13 CO+, we fit the 
DCN/HCN ratio over the whole disk and determine the 
DCN/HCN ratio to be 1.7±0.5 x 10~ 2 . To again empha- 
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Fig. 6. — Iso-x 2 surfaces of (Rout,Pi) for HCO+, DCO+ (as in 
Model 2) and HCN. Contours correspond to the 1 to 6 a errors. For 
DCO" 1 ", the index values larger than 3 at around 90 AU indicate 
the ratios of DCO+/HCO+ larger than 1, so the \ 2 surface is not 
calculated beyond that. 

size the impact of the assumed vertical distribution on 
the derived column densities, the DCN/HCN ratio could 
be as high as 5 x 10~ 2 if HCN and DCN are distributed 
vertically over the same region as is HCO + . 

Highly fractionated DCN/HCN ratios have been mea- 
sured in comets. In the coma of comet Halc-Bopp, for ex- 
ample, Meier et al. (1998) reported the ratio to be around 
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Fig. 7. — The beam-averaged H 13 CO + J=4-3 spectra at the 
continuum (ste liar) position. The SMA data are presented by the 
solid histogram, and the simulated model by the dashed histogram. 
The vertical dotted line indicates the position of the fitted Vlsr 
for HCO+ J=3-2. 
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Fig. 8. — Models of DCO+ with different distributions of radial 
column densities. 

2.3 x 1CT 3 , but higher DCN/HCN ratios - (D/H) H CNjet 
«0.025 are detected from the pristine material sublimed 
from icy grains ejected in jets from the nucleus which 
may present a more representative sampling of cometary 
ices that have not experienced significant thermal pro- 
cessing (Blake et al. 1999). Such ratios are consistent 
with those found here in the TW Hydrae disk, indicat- 
ing that high D/H ratios in comets could originate from 
material in the outer regions of disks where in situ deu- 
terium fractionation is ongoing, rather than requiring an 
inheritance from interstellar material. 

4.3. Upper limits for H 2 D+ , D 2 H+ and EDO 

In the disk midplane, H2 is expected to be gaseous 
and the molecular ion formed by the cosmic ray ioniza- 
tion of H2, H^, is expected to be the most abundant 
ion. Unfortunately, H^" is only detectable in cold gas 
via infrared absorption. In its deuterated forms, how- 
ever, H2D + and even D2H + are expected to be abun- 
dant in the cold, dense gas (Ceccarelli & Dominik 2005). 
The ground-state transition of ortho-H2D + was first de- 
tected in a young stellar object (NGC 1333 IRAS 4A) 
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Fig. 9. — DCO + J=3-2 channel maps toward TW Hydrae and the simulated model distributions depicted in Figure 8. 



by Stark et al. (1999) and in a prestcllar core (L1544) 
by Caselli et al. (2003). Both H 2 D+ and D 2 H+ have 
been detected toward another pre-stcllar core 16293E via 
their ground-state submm rotational lines (Vastel et al. 
2004) . The inclusion of multiply deuterated H J in chem- 
ical models leads to predictions of higher values of the 
D/H ratio in cold, high-density regions of the interstellar 
medium. Similarly, in the dense, cold disk midplane, CO 
is depleted, and high abundances of H 2 D + and D 2 H + 
are expected. For this reason, Ceccarelli et al. (2004) 
searched for the ground-state transition of ortho-H 2 D+ 
with the Caltech Submillimeter Observatory (CSO), and 
reported 3.2er and 4.7a detections toward TW Hydrae 
and DM Tauri, respectively. With the 400 GHz and 690 
GHz receiver-equipped SMA antennas, we searched for 
the 372 GHz ortho-H 2 D+ l M -li, and 690 GHz para- 
D 2 H + lio-loi lines toward TW Hydrae. No significant 
emission signals were detected. Here we discuss the up- 
per limits and their implications. 

Our 3-antenna SMA observations give a la upper 
limit for the ortho-H 2 D + li,i-li,o line emission of 1.2 
Jy beam -1 km s -1 with a 4'.'7 x 3'.'8 synthesized beam. 
To compare the result with single dish data, the extent 
of the source emission must be known. Since H 2 D + was 
observed along with the H 13 CO + 4-3 line in a dual- 
receiver observation on 28 December 2006 and H 13 CO + 



4 3 has been clearly detected at JCMT (van Dishoeck 
et al. 2003), we can compare the intensities of these two 
lines between the SMA data and the single dish observa- 
tions to constrain the emitting areas. The H 13 CO + 4-3 
line was detected at JCMT with an integrated intensity 
of 0.07 K km s -1 in a 13" beam (van Dishoeck et al. 
2003); while for the SMA the integrated intensity of this 
line is determined to be 0.61 Jy beam -1 km s -1 with 
a beam of 4'.'1 x l'/8. If the extent of the emission of 
H 13 CO + is similar to that of H 2 D + , our interferometric 
upper limit of Jy beam -1 km s -1 for H 2 D + (consider- 
ing the small change of the beam sizes) becomes 0.09 
K km s -1 (la) or 0.27 K km s -1 (3 a) upper limits, 
which is consistent with 2a upper limits of about 0.2 
K km s -1 by the JCMT (Thi et al. 2004) but slightly 
lower than the 3.2a detection of 0.39 K km s -1 by Cec- 
carelli et al. (2004). We estimate the la upper limit of 
the ortho-H 2 D+ column density to be 1.7 x 10 12 cm -2 
according to the Equation 4 of Vastel et al. (2004), which 
is slightly less than the 2a upper limit estimate of 4.4 x 
10 12 cm -2 by Thi et al. (2004) since the SMA data had 
a smaller noise level. Of course this analysis assumes the 
extent of H 2 D+ is similar to that of H 13 CO+. With the 
deployment of more 400 GHz receivers on the SMA and 
further H 2 D + observations, it should be possible to pro- 
vide rather better constraints on the H 2 D + abundance 
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Fig. 10.— Top: Velocity channel map of the DCO+ J=3-2 toward TW Hydrae. The angular resolution is 2'.'6 X T.'6 at PA 2.8°. The 
cross indicates the continuum (stellar) position. The axes are offsets from the pointing center in arcseconds. The la contour step is 0.12 Jy 
beam -1 and the contours start at 2a. Middle: channel map of Model 3 with the same contour levels. Bottom: difference between Model 
3 and data on the same contour scale. 
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Fig. 11.— Top: Velocity channel map of the HCN J=3-2 toward TW Hydrae. The angular resolution is 1'.'6 X l'.'l at PA -0.5°. The 
cross indicates the continuum (stellar) position. The axes are offsets from the pointing center in arcseconds. The la contour step is 0.35 
Jy beam -1 and the contours start at 2a. Middle: channel map of the best-fit model with the same contour levels. Bottom: difference 
between the best-fit model and data on the same contour scale. 
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Fig. 12. — Top left: DCN J=3-2 velocity channel maps (red: 2.84 km s_1 > blue: 2.28 km s -1 ) from TW Hydrae, overlaid on the 217 
GHz dust continuum map (gray scale). The cross indicates the position of the continuum peak. Top right: the simulated model for DCN 
J=3— 2. The la contour step is 0.09 Jy beam -1 and the contours start with la. Bottom: the beam-averaged DCN 3—2 spectra at the 
continuum (stellar) position. The SMA data are presented by the solid histogram, and the simulated model by the dashed histogram. The 
vertical dotted line indicates the position of the fitted Vlsr f° r HCN J=3— 2. 
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in the disk. 

We estimate the la upper limit for para-D 2 H + lio loi 
to be 5.35 Jy beam -1 km s _1 with a beam of 3'.'3 x l'/3. 
The la upper limit to the para-D 2 H+ column density is 
estimated to be 9.0 x 10 14 cm~ 2 . This is less constrained 
than H 2 D+ due to the relatively poor system sensitivity 
at 690 GHz. 

For HDO 3i,2-2 2 ,i, the la upper limit is 0.10 
Jy beam -1 km s _1 . Assuming an excitation temperature 
of 30 K, the upper limit for the HDO column density is 
2.0 x 10 14 cm~ 2 . Since the lower state energy level of 
this line is nearly 160 K, it must originate from warm 
regions of the disk which are quite distinct from the cold 
layers where HDO ground state transition absorption, as 
found in DM Tauri (Ceccarelli et al. 2005), must arise - 
although we note that the detection of the HDO absorp- 
tion line in DM Tauri has been disputed by Guillotcau 
ct al. (2006). 

5. SUMMARY AND CONCLUSIONS 

Observations of dcuterated species in circumstellar 
disks are important to understand the origin of primitive 
solar system bodies in that they can directly constrain 
the deuterium fractionation in the outer regions where 
cometary ices are likely formed. Spatially resolved obser- 
vations of the D/H ratios in disks enable the comparison 
of the fractionation measured in comets such as Hale- 
Bopp (Blake et al. 1999) with the specific conditions at 
each disk radius. We have presented the first images 
of DCO + and DCN emission from the disk around a 
classical T Tauri star, TW Hydrae, along with images 
of the HCN and HCO+ J=3-2 lines. The observations 
of deuterium fractionation serve as a clear measure of 
the importance of low-temperature gas-phase deuterium 
fractionation processes. These observations strongly sup- 
port the proposed link among high gas densities, cold 
temperature and enhanced deuterium fractionation. Dc- 
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tailed chemical models are still needed to explain how 
DCO + disappears from the outer part of the disk. 

The similarity of the D/H ratios in cold clouds, disks 
and pristine cometary material has been used to argue 
that the gas spends most of its lifetime at low tempera- 
tures and is incorporated into the disks before the enve- 
lope is heated, i.e. before the Class I stage. By combining 
self-consistent physical models and 2D radiative transfer 
codes to interpret high spatial resolution millimeter-wave 
molecular images, we are only now beginning to investi- 
gate the radial and vertical distributions of molecules in 
disks. The radial distribution of DCO + in the disk of 
TW Hydrae indicates that in situ deuterium fractiona- 
tion is ongoing. The molecular evolution within disks 
must therefore be considered in the investigation of the 
origin of primitive solar system bodies. 

We have obtained less stringent constraints on the ver- 
tical distributions of molecules in the disk of TW Hydrae. 
To address the ambiguity present in the analysis of single 
objects, data from a robust sample of disks is needed, in 
particular one that covers a range of disk inclinations. 
More sensitive observations are also needed for the rare 
isotopologues H 13 CN, H 13 CO+ and, of course, DCN, to 
understand the radial and vertical gradient of deutcratcd 
species in these disks. In the future, observations of DCN 
and other species with the Atacama Large Millimeter Ar- 
ray will provide further insight into the chemical state of 
protoplanetary disks. 
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